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ABSTRACT 

Context. HD 181327 is a young main sequence F5/F6 V star belonging to the j 3 Pictoris moving group (age ~ 12 Myr). It harbors an optically thin 
belt of circumstellar material at radius ~90 AU, presumed to result from collisions in a population of unseen planetesimals. 

Aims. We aim to study the dust properties in the belt in details, and to constrain the gas-to-dust ratio. 

Methods. We obtained far-infrared photometric observations of HD 181327 with the PACS instrument onboard the Herschel Space Observatory* , 
complemented by new 3.2 mm observations carried with the ATCA** array. The geometry of the belt is constrained with newly reduced 
HST/NICMOS scattered light images that allow the degeneracy between the disk geometry and the dust properties to be broken. We then use 
the radiative transfer code GRaTeR to compute a large grid of models, and we identify the grain models that best reproduce the spectral energy 
distribution (SED) through a Bayesian analysis. We attempt to detect the oxygen and ionized carbon fine-structure lines with Herschel/PACS spe- 
ctroscopy, providing observables to our photochemical code ProDiMo. 

Results. The HST observations confirm that the dust is confined in a narrow belt. The continuum is detected with Herschel/PACS completing 
nicely the SED in the far-infrared. The disk is marginally resolved with both PACS and ATCA. A medium integration of the gas spectral lines only 
provides upper limits on the [Oil and [CII] line fluxes. We show that the HD 181327 dust disk consists of micron-sized grains of porous amorphous 
silicates and carbonaceous material surrounded by an important layer of ice, for a total dust mass of ~0.05 M e (in grains up to 1 mm). We discuss 
evidences that the grains consists of fluffy aggregates. The upper limits on the gas atomic lines do not provide unambiguous constraints: only if 
the PAH abundance is high, the gas mass must be lower than ~17 M e . 

Conclusions. Despite the weak constraints on the gas disk, the age of HD 181327 and the properties of the dust disk suggest that it has passed the 
stage of gaseous planets formation. The dust reveals a population of icy planetesimals, similar to the primitive Edgeworth-Kuiper belt, that may 
be a source for the future delivery of water and volatiles onto forming terrestrial planets. 

Key words, stars: individual: HD 181327 - circumstellar matter - infrared: planetary systems - radiative transfer 


1. Introduction 

The evolution of planets in the solar system is intimately con- 
nected to the existence of a reservoir of planetesimals in its outer 
regions. The depletion of the Kuiper belt and the main aster- 
oid belt that occurred during the late heavy bombardment not 
only delivered large amounts of water, volatiles and carbona- 
ceous material onto the inner planets, but it also reduced the 
later rate of catastrophic impacts onto the Earth, opening the 
path to the emergence of life (Morbidelli et al. 2000; Greaves 
& Wyatt 2010; Hartogh et al. 2011). Although detecting Kuiper 
belt-like objects around nearby planetary systems will remain 
an unachievable goal for the foreseeable future, their colli- 
sional erosion produces circumstellar disks of dust responsible 
for a characteristic excess emission detectable at infrared and 
sub-millimeter wavelengths. Recent studies using the Herschel 
Space Observatory reveal that as much as ~30% of nearby Main- 
Sequence F, G, K stars are surrounded by cold debris rings anal- 
ogous to the Kuiper belt (latest Herschel / DUNES open time 
key program results, Eiroa et al., in prep.). However, the cold 

* Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and with im- 
portant participation from NASA. 

** The Australian Telescope Compact Array is operated by the 
Australian Telescope National Facility (ATNF) managed by CSIRO. 


temperature and the large grain sizes make it impossible to de- 
tect solid-state features in their spectrum, thus preventing unam- 
biguous identification of the dust composition. We are thus left 
with model-dependent methods to constrain the dust composi- 
tion from color and polarimetric measurements, and from the 
spectral energy distribution (hereafter SED). 

Another uncertainty in planet formation theories lies in the 
lack of detailed observational constraints for the mutual dust and 
gas dissipation timescales in circumstellar disks, or more pre- 
cisely, on how the gas-to-dust mass ratio evolves with time and 
location in disks at the very early stages of planet formation. The 
study of stellar clusters of different ages shows that, statistically, 
a population of young stars loses its massive dust disks in only 
a few million years on route to the main sequence (Carpenter 
et al. 2005), although, for individual objects, the transition from 
an optically thick to an optically thin disk is expected to occur 
more quickly (a million years or less, Cieza et al. 2007; Currie 
& Sicilia- Aguilar 201 1). However not much is known about the 
characteristic evolutionary timescales for the gas in these disks. 
Observing young stars at different evolutionary stages can help 
elucidate the extent to which gas and dust simultaneously dissi- 
pate in disks. This requires observations of both the dust and gas 
components and, at the same time, modeling of the continuum 
and line emission. 


Article published by EDP Sciences 
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Table 1. Star and disk properties. 


Parameter 

Value 

Reference 

Spectral type 

F5/F6V 

Nordstrom et al. (2004) 

Magnitude V 

7.0 

Bessell (2000) 

Gravity (log(g)) 

4.510 

Chen et al. (2006) 

Distance (r/*) 

51.8 pc 

Holmberg et al. (2009) 

Age 

12 Myr 

Zuckerman & Song (2004) 
Mamajek et al. (2004) 

Luminosity (L*) 

3.33 L® 

this study 

Mass (M*) 

1.36 M 0 

this study 

Disk PA (b) 

107 ± 2° 

Schneider et al. (2006) 

Disk inclination 

31.7 ± 1.6° 

Schneider et al. (2006) 

7-ir /L* 

~2 x 10“ 3 

this study 


Notes. (a) From integrating the synthetic stellar spectrum. (i,) Position 
angle, east of north. 


As part of the GASPS (Gas in Protoplanetary Systems) 
Herschel Open Time Key Programme, we observed HD 181327, 
a young (~12Myr) F5.5V star located at 51.8pc (Holmberg 
et al. 2009, Table 1) belonging to the p Pictoris moving group 
(Zuckerman & Song 2004; Mamajek et al. 2004). HD 181327 
was identified as a debris disk hosting main sequence star 
through its SED with a fractional infrared luminosity Tir/L* = 
0.2% (Mannings & Barlow 1998). Schneider et al. (2006), from 
NICMOS coronagraphic observations in scattered light, dis- 
covered a ring-like disk of circumstellar dust located at radius 
~89 AU from its star. The ring is inclined by 31.7 ± 1.6° from 
face-on and shows an apparent brightness assymetry with re- 
spect to the minor axis that is well explained by a strong di- 
rectionally preferential scattering. No photocenter nor pericenter 
offset is seen in the ring relative to the position of the central star 
and the disk appears axisymmetric. Recent reprocessing of these 
data (Schneider et al., in prep.) confirms the annulus shape of 
the disk, and suggests a narrow radial distribution of the dust in 
the system (FWHM = 24.5 AU). A low surface brightness dif- 
fuse halo is seen in the NICMOS image and in complementary 
HST/ACS images at a distance of ~4" < r < 9" from the cen- 
tral star, which may correspond to a population of very weakly 
bound grains originating in the main ring. Mid-infrared images 
at 18.3 gm by Chen et al. (2008) suggest azimuthal asymmetries 
in the density profile of the ring with respect to the apparent mi- 
nor axis. These 18.3/rm observations also show that little to no 
emission (<15%) comes from inside the ring resolved in scat- 
tered light. Therefore, the density profile of the HD 181327 dust 
disk should be well constrained by the scattered light and ther- 
mal emission images, allowing a detailed study of the dust prop- 
erties. Interestingly, from IRS and MIPS-SED data, Chen et al. 
(2008) also inferred the presence of crystalline ice in the debris 
disk. 

In this paper we perform SED and line fitting of new 
HerschelfPACS observations of HD 181327 (presented in 
Sect. 2), supplemented by data from the literature. We also 
present new ATCA imaging of the disk (Sect. 2). In Sect. 3, 
we derive surface density profiles from the newly reduced 
HST/NICMOS data, which well constrain the ring geometry. We 
describe our dust model and present the results of the SED fitting 
in Sects. 3 and 4. We assess the gas content in this young debris 
disk in Sect. 5. Finally we discuss the results in Sects. 6-8 and 
present our conclusions in Sect. 9. 


Table 2. Photometry results for HD 181327 from chop-nod 1 and scan 
map 2 observations. 


Wavelength 

(pm) 

/A 

* cont 

(Jy) 

Absolute flux cal. 
uncertainty 

70 1 

1.827 ±0.0069 

10% 

100 2 

1.337 + 0.0082 

10% 

160 1 

0.767 ± 0.015 

20% 

160 2 

0.772 + 0.011 

20% 


Notes. (tl The ±lcr error bars on the flux measurements include statis- 
tical uncertainty only. 


2. Observations and data reduction 

2.1. Herschel/PACS continuum 

We obtained new far-IR photometry and spectroscopy of 
HD 181327 using the PACS instrument onboard the Herschel 
Space Observatory (Pilbratt et al. 2010; Poglitsch et al. 2010). 
The observations were carried out in two modes: 1) 159 s of 
point-source chop-nod mode imaging at 70 and 160 gm (obsid: 
1342183658) and 2) two 276 s scan map imaging observations 
at 100 and 160 /rm (obsids: 1342209057-8). The scan maps were 
executed with the medium scan speed (20" /s), at scan angles of 
70° and 1 10°. The two cross scan maps at different scan angles 
were combined to increase the signal-to-noise by cutting down 
on streaking in the image background along the scan direction. 

We reduced the photometry with the Herschel Interactive 
Processing Environment (HIPE, Ott 2010) version 4.2. The un- 
certainties in the absolute flux calibration for this version are 
10% at 70 and 100 /jm and 20% at 160 /jm (see the PACS 
Scan Map Astronomical Observation Templates release note 1 ). 
For the 70" and 160 pm chop-node images, we performed aper- 
ture photometry using apertures of 7" and 1 1", respectively. For 
the 100 and 160 pm scan map observations, we used apertures 
of radii 12" and 24", respectively. An annulus for sky back- 
ground estimation was placed from 10"-20" beyond the aper- 
ture in each case, and an aperture correction was applied to the 
measured fluxes. Both estimations at 160 /jm are fully compat- 
ible within the error bars; the 160 pm scan-map flux is used in 
the rest of the study. Since the source is marginally resolved at 
70 and 100 pm (see Sect. 2.1), we calculated our own aperture 
corrections following the procedure described in Poglitsch et al. 
(2010). The photometry results are presented in Table 2. 

HD 1 8 1 327 is more extended than the Herschel point-spread 
function (PSF) in both the 70 gm and 100 gm images. The 70 pm 
image of HD 181327 appears in Fig. 1, along with the 70 p m 
image of a PSF reference star (HD 148387) observed using the 
same chop-nod mode. The upper panel of Fig. 2 shows the 
azimuthally-averaged radial brightness profiles for HD 181327 
and HD 148387. In the case of the 100 gm scan map observa- 
tion of HD 181327, we used a different PSF reference observa- 
tion, i.e. a scan map of a Boo. The 100 pm a Boo radial profile 
appears in Fig. 2, bottom panel. At both wavelengths, the az- 
imuthally averaged radial profile of HD 181327 is clearly ex- 
tended beyond that of the PSF reference. The FWHM of the 
HD 181327 70 gm profile is 6744, compared to a FWHM of 
5728 for HD 148387. At 100 ^m, the FWHM for HD 181327 is 
77 5 8, compared to 6794 for a Boo. 


1 Available at 

http: //herschel . esac . esa. int/AOTsReleaseStatus . shtml 
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Fig. 1. PACS 70 pm imaging of {left) HD 181327 and {right) a PSF ref- 
erence star, HD 148387. Both images were acquired in chop-nod mode. 
The field of view of each image is 25" x 25". Brightness contours at 5, 
10, 25, 50 and 75 times the background noise in the HD 181327 image 
are overlaid on both images. 



Radius (arcsec) 


Fig. 2. Radial profiles of HD 181327 at 70 pm {top), compared to 
HD 148387 as a PSF reference, and at 100 pm {bottom), compared to 
a Boo as a PSF reference. Both plots show the azimuthally averaged, 
normalized radial brightness profile vs. radius from the star. Each error 
bar is the standard deviation of the brightness values at that radius. The 
FWHM values for the profiles are given on the plot. 


2.2. ATCA 3.2 mm observations 

The Australia Telescope Compact Array (ATCA) was used to 
observe HD 181327 at 3.2 mm in May and September 2010. 
The observations were carried out in continuum double sideband 
mode, with 2046 x 1 MHz channels per sideband with an ef- 
fective total bandwidth of just over 2 GHz. The ATCA has five 
antenna with 3-mm receivers, and our observations were carried 
out with the compact H75 array configuration with baselines of 
31-98 m. The quasars PKS 1815-553 and PKS 1253-055 were 
used for the gain and bandpass calibrations respectively, and 
Uranus was used for the absolute flux calibration. The ATCA 
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RA (arcsec) 


Fig. 3. ATCA 93.969 GHz (3.192 mm) map of HD 181327 with con- 
tours at -2 (dashed red line), 2 and 3cr (solid blue lines), where cr is the 
image rms given by 1.1 x 10 -4 Jy/beam. The dashed blue area repre- 
sents the beam. The Gaussian fit is overplot with contours at 2 and 3cr 
(solid green lines). The red and green crosses denote respectively the 
star position and the Gaussian fit center. 


calibration uncertainty is estimated to be about 20% in the 3-mm 
band. The data were calibrated and reduced with the MIRIAD 
package (Sault et al. 1995). The frequency pair for the obser- 
vations were centered at 93 and 95 GHz, and during the re- 
duction process the data was averaged over 64 channels to im- 
prove the signal-to-noise, shifting the frequency pairs to 92.968 
and 94.969 GHz. The final CLEANed map combined the two 
sidebands, resulting in an image at 93.969 GHz (3192.5 pm) 
presented in Fig. 3. Our total on-source integration time was 
770 minutes. We determined an integrated flux of 7.2 x 10 -4 Jy, 
with an rms of cr = 1.1 x 10 -4 Jy/beam. Including the 20% cal- 
ibration uncertainty, this yields an uncertainty of 2.5 x 10 -4 Jy. 
The disk is marginally resolved. Assuming a Gaussian pro- 
file, we find a PA of ~102°, consistent with that seen on the 
NICMOS image (107±2°, Table 1), and along which the FWHM 
is ~7.8" (beam: 3.5" x 2.8", PA = 85.4°). The center of the 
Gaussian fit is also shifted by 0.6" toward the West with re- 
spect to the star position 2 . Structures consistent with a belt at 
~100 AU are detected along the semi-minor axis with a peak 
flux of 4.01 x 10 -4 Jy/beam on the NW side. We note a possible 
brightness asymmetry with respect to the SE side. 

2.3. Herschei/PACS spectroscopy 

We obtained far-infrared spectroscopy using the PACS integral 
field unit which provides a 5 x 5 array of 9'.' 4x 9'.'4 spaxels. The 
data were taken in chop-nod mode using the line scan and range 
scan modes. The spectra cover eight small wavelength ranges 
(~1.8 pm wide) targeting particular emission lines, including 
[Oi] at 63 pm and [CII] at 158 pm. The spectral resolution of 
the data varies between R ~ 1 100 and 3400. Additional details 
on the observations appear in Mathews et al. (2010). 

The data were calibrated with the Herschel Interactive 
Processing Environment (HIPE) version 4.2. The wavelength- 
dependent flux and aperture corrections provided by the PACS 
development team were applied, giving an uncertainty in the 

2 The pointing accuracy of ATCA is ~ 1". 
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Fig. 4. PACS spectrum of HD 18 1327, centered on the wavelength of 
the [O i] 63 pm emission line. The linear fit to the continuum flux is 
overplotted with the red line. The [Ot] 63 pm emission line was not 
detected. 


Table 3. Spectroscopy results for HD 181327. 


Line 

Wavelength 

(pm) 

cont 

(Jy) 

F* 

int 

(xlO -18 W/m 2 ) 

[Ot] 

63.1852 

1.731 ± 0.236 

<9.82 

CO 36-35 

72.8429 

1.491 ±0.185 

<18.41 

o-H 2 0 

78.7414 

1.431 ±0.187 

<10.51 

CO 29-28 

90.1630 

1.264 ±0.222 

<8.90 

[Ot] 

145.5350 

0.804 + 0.107 

<8.48 

[Cn] 

157.6790 

0.767 + 0.175 

<7.02 

o-H 2 0 

179.5265 

0.388 + 0.201 

<6.87 

DCO + 22-21 

189.5700 

0.581 + 0.378 

<16.71 


Notes. m The ±l<x error bars on the flux measurements include statis- 
tical uncertainty only. The absolute calibration uncertainty amounts to 
30%. ® 3cr upper limit. 


absolute flux calibration of 30% (see the PACS Spectroscopy 
Performance and Calibration Document 2 ). Significant emission 
was detected only in the central spaxel. No emission lines are 
seen, only continuum flux. The central spaxel spectrum centered 
on the [O i] 63 pm line appears in Fig. 4. 

The continua were binned at half the instrumental FWHM 
and fitted with 1st degree polynomials, using the fluxes and er- 
rors produced by our HIPE reduction script. The HIPE uncer- 
tainties appeared to be slightly too small, however. New flux 
random errors were calculated from the rms of the data minus 
the continuum fits, over wavelength ranges spanning 6 to 14 
times the instrumental FWHM, centered on the expected line 
rest wavelengths. The continuum fluxes at the rest wavelengths 
of the expected emission lines appear in Table 3. To determine 
the upper limit on the integrated line flux (Fj nt ) for each avail- 
able emission line, we integrated the data over a small wave- 
length range, either ±1 or ± 1 .5 resolution elements, centered on 
the rest wavelength of the expected emission line. The flux errors 
were propagated to determine the l<x upper limit on each inte- 
grated line flux. The 3cr integrated line flux upper limits appear 
in Table 3 as well. 



Distance from the star in All 

Fig. 5. New HST/NICMOS 1.1 /im surface brightness profiles (top, 
nominal and ±lcr profiles), and surface density profiles resulting from 
direct inversion of the brightness profiles (bottom), assuming two pos- 
sible values for the asymmetry parameter. The color-coding is the same 
for both panels. The striped areas identify regions dominated by PSF- 
substraction residuals where the density has been extrapolated. 


3. Model setup 

The dust and gas models developed in Sects. 4 and 5 rely on 
some assumptions on the spatial distribution of the dust and the 
gas, as well as on the dust properties. These are detailed below in 
Sects. 3.1 and 3.4. Additional data collected from the literature 
and stellar properties are summarized in Sect. 3.2. 

3. 1. New surface density profile for the HD 181327 ring 

Our spatially resolved imaging constrains the location of the 
dust. We advantageously impose this constraint in order to ex- 
amine the dust properties in details and to get an upper limit 
on the gas mass from the Herschel non-detections of gaseous 
emission lines. The disk around HD 181327 was spatially re- 
solved in scattered light with HST/NICMOS 1 . 1 pm imaging by 
Schneider et al. (2006), which alleviates many uncertainties on 
the disk structure. The NICMOS images revealed a ring-like disk 
with a peak density position ro = 1.71” (~88. 6 AU, according to 
the new distance to the star, 51.8 pc, by Holmberg et al. 2009), 
with an azimuthally medianed photometric FWHM of ~0.538" 
(~28 AU). 

Recently, Schneider et al. (in prep.) improved the data reduc- 
tion of the original HST/NICMOS observations of HD 181327 
thanks to PSF subtraction template images better matched to the 
stellar PSF structure in the HD 181327 images. This results in 
an overall improvement of the signal-to-noise, especially inside 
the ring, and leads to a slight sharpening of the belt. The new 
profiles, displayed in the top panel of Fig. 5, show a peak at a 
distance of 88.6 AU and a surface brightness fainter by ~21% 
(0.76 mjy arcsec -2 at the peak position) compared to the origi- 
nal Schneider et al. (2006) study. 

We derive the surface density profile Z(r) by inverting the 
new scattered light surface brightness profile. We employ the 
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methodology described in Augereau & Beust (2006), assum- 
ing a disk inclination of i = 31.7° from face-on. The inversion 
procedure yields the product <x sca Z(r), where <x sca is the mean 
scattering cross section of the grains at 1 . 1 gm. This is equiva- 
lent to the product of the vertical optical depth T±(r), multiplied 
by the mean albedo u> at that wavelength (cr sca Z(r) = u>T ± (r), 
see Eqs. (2), (3) in Schneider et al. 2006). The inversion solu- 
tion depends on both the anisotropic scattering phase function, 
and on the disk vertical profile. Nevertheless, Schneider et al. 
(2006) showed that the density profile is largely independent of 
the shape for the vertical dust distribution in the range of ac- 
ceptable vertical scale heights II g. We here assume a Gaussian 
vertical profile with a scale height Hg = 6 AU at the observed 
peak radius of the ring rg = 88.6 AU {Hg/rg ~ 0.068, consistent 
with the theoretical lower limit of 0.04 from Thebault 2009). The 
observed phase function of the disk at 1 . 1 pm is well-represented 
by a Henyey-Greenstein scattering phase function with asymme- 
try parameter I^hgI = 0.30 ± 0.03. This value was discussed by 
Chen et al. (2008) in light of asymmetries identified in T-ReCS 
thermal emission images in the g fl -band. Chen et al. (2008) sug- 
gest that the asymmetries may have led to an overestimate of 
|{7 hgI, which may be as small as 0.16 ± 0.04. We thus perform 
the inversion for both |gr H d = 0-16 and |gr H Gl = 0.30. 

The results are displayed in the bottom panel of Fig. 5, where 
we also show the impact of the surface brightness uncertainties 
on the inferred surface density profile. We find an actual surface 
density peak at 89.5 AU from the star. A higher I^hgI leads to 
a globally higher surface density profile. The surface brightness 
measured interior to ~1.5" (~78 AU) is not reliable because of 
PSF subtraction residuals. We fit the profiles from 79 to 83 AU 
obtained for both I^hgI values by a very steep power-law func- 
tion of the form ~r +10±l . We extrapolate this profile interior to 
78 AU as shown in Fig. 5 (striped area). This extrapolation is 
supported by the mid-IR results of Chen et al. (2008), who show 
that the inner part of the disk is largely cleared toward the center, 
containing little to no dust. 

Outside of the belt peak position, the improved data reduc- 
tion results in a profile that falls off as ~r ~ 4J±03 from ~92 
to 162 AU, a steeper fall-off than the Schneider et al. (2006) 
original profiles (~r -3 ). The dust is thus more tightly confined 
to the parent-belt than previously estimated, and the profile is 
also steeper than theoretical expectations. The exact value of 
the outer slope can only slightly impact the modelling results, 
as most of the emission arises from the inner part of the belt 
where the temperature is larger, except possibly for the longer- 
wavelengths. 

Eventually, the differences between the several radial profiles 
in Fig. 5 essentially translate into a scaling factor that impacts the 
total disk mass estimate. We estimate that the uncertainty on the 
asymmetry parameter adds an additional ~4% to the uncertainty 
on our final mass values. The density profile obtained for |<7 H gI = 
0.30 is used in the rest of this study. We stress that, in order to 
focus on the other disk properties, we make the assumption that 
this profile is valid for any dust grain constrained by the SED, i.e. 
for grain sizes ranging from the micron to the millimetre scale. 
This assumption assumes no radial segregation in dust sizes and 
it remains reasonable as long as the flux at any wavelength arises 
predominantly from regions close to the density peak position 
(see Sect. 6.1). 

3.2. Spectral energy distribution 

The Herschel and ATCA data presented in this paper are com- 
plemented by IR and sub-mrn data collected from the literature. 


Table 4. HD 181327 photometric observations from the literature. 


Wavelength 

(urn) 

Flux 

(Jy) 

Uncertainty 

(Jy) 

Instrument 
& Reference 

12 

0.114 

0.011 

IRAS [1] 

15.0 

0.094 

0.004 

Spitzer/ IRS [21 

16.3 

0.097 

0.005 

Spitzer/ IRS [2] 

17.4 

0.107 

0.005 

Spitzer/IRS [21 

18.3 

0.114 

0.003 

T-ReCS [3] 

18.7 

0.125 

0.007 

Spitzer/ IRS [21 

19.8 

0.144 

0.009 

Spitzer/ IRS [21 

20.4 

0.149 

0.008 

Spitzer/ IRS [21 

21.1 

0.157 

0.020 

Spitzer/IRS [21 

23.1 

0.224 

0.007 

Spitzer/ IRS [21 

23.7 

0.223 

0.009 

Spitzer /mPS [11 

25 

0.306 

0.024 

IRAS [11 

25.7 

0.303 

0.014 

Spitzer/IRS [2[ 

28.0 

0.384 

0.014 

Spitzer/ IRS [2[ 

30.6 

0.488 

0.026 

Spitzer/IRS [2[ 

33.1 

0.623 

0.025 

Spitzer/IRS [2[ 

35.5 

0.718 

0.034 

Spitzer/IRS [2[ 

54.5 

1.88 

0.09 

MIPS-SED [41 

58.0 

1.83 

0.07 

MIPS-SED [41 

60 

1.73 

0.17 

ISO [51 

61.4 

1.98 

0.08 

MIPS-SED [41 

64.9 

2.01 

0.08 

MIPS-SED [41 

68.4 

1.98 

0.08 

MIPS-SED [41 

71.4 

1.73 

0.12 

Spitzer/ MIPS [1] 

71.6 

1.99 

0.09 

MIPS-SED [41 

75.1 

1.93 

0.10 

MIPS-SED [41 

78.4 

1.81 

0.10 

MIPS-SED [41 

81.8 

1.72 

0.11 

MIPS-SED [41 

85.3 

1.62 

0.11 

MIPS-SED [41 

88.8 

1.47 

0.11 

MIPS-SED [41 

90 

1.41 

0.14 

ISO [51 

92.0 

1.46 

0.13 

MIPS-SED [41 

95.5 

1.37 

0.17 

MIPS-SED [41 

100 

1.7 

0.2 

IRAS [1] 

155.9 

0.77 

0.09 

Spitzer /MARS [ 1 1 

170 

0.736 

0.192 

ISO [51 

870 

0.052 

0.006 

LABOCA [6] 

3190 

7.2 x 10“ 4 

2.5 x 10“ 4 

ATCA [7] 


Notes. The Herschel continuum values used for the dust modeling are in 
Tables 2 and 3. All uncertainties are lcr and include absolute calibration 
error. 

References. [11 Schneider et al. (2006), [21 this study, data reduced with 
the Spitzer/clA team pipeline (Lahuis et al. 2006), [3] Gemini, Chen 
et al. (2008), [4] Spitzer, Chen et al. (2006), [51 Moor et al. (2006), [6] 
APEX, Nilsson et al. (2009), [7] this study, Sect. 2.2. 


These are listed in Table 4. The shortest wavelength flux we use 
is the 12 pm IRAS photometry for which the emission shows 
no excess from the stellar photosphere level. This observation is 
documented in Schneider et al. (2006) together with IR-excess 
from the IRAS Faint Source Catalogue and obtained with the 
Spitzer / MIPS instrument. These are complemented by ground- 
based mid-IR measurements performed with the T-ReCS imager 
on Gemini-S by Chen et al. (2008), and 13 points extracted from 
the 5-35 pm Spitzer/IRS spectrum that we reduced using the c2d 
Fegacy team pipeline (Lahuis et al. 2006). 

Additional MIPS observations were obtained in SED-mode 
by Chen et al. (2006) between 54.5 and 95.5 pm, well- 
constraining the shape of the disk SED around its maximum 
emission wavelength. Far-IR excess was reported by Moor et al. 
(2006) using the ISO database. The shape of the SED in the 
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Fig. 6. Observed spectral energy distribution of HD 181327 (see 
Tables 2-4). 


far-IR domain (100-200 pm) is well constrained by our new 
Herschel measurements. 

Because the uncertainties on the dust optical constants at 
millimetre wavelengths are high, we chose not to use the ATCA 
flux to constrain the dust model (see Sect. 6.1 for a discus- 
sion on the compatibility of the model to these observations). 
Therefore, the sub-mrn regime is only represented by a measure- 
ment with the LABOCA bolometer array at the 12-m telescope 
APEX (Nilsson et al. 2009, A ~ 870 pm), which thus represents 
a key piece of information. These archival data, together with 
our new observations, constitute the 49 measurements that are 
displayed in Fig. 6. 

3.3. Star properties and dik geometry 

We derive excesses from a comparison of the SED with a syn- 
thetic stellar spectrum extracted from the “NextGen” model at- 
mosphere grid (Hauschildt et al. 1999) and scaled to match 
the visible magnitude V = 7.0 mag (Bessell 2000), assuming 
log(gr) = 4.5, T e ff = 6600 K, which leads to a stellar luminosity 
of 3.33 Lq and a mass of 1 .36 M Q . This stellar model shows good 
agreement with the Spitzer / IRS spectrum at wavelengths shorter 
than 10-15pm, longward of which the excess starts dominating 
over the stellar photospheric emission. The star properties and 
disk geometry adopted in this study are summarized in Table 1 . 

3.4. Dust properties: size distribution and composition 

For both the dust and gas modeling, we make the assumption 
that the grain population has a single differential size distribu- 
tion valid anywhere in the ring. This is a central hypothesis that 
necessarily limits the uniqueness of the solution. We discuss its 
impact in Sect. 6. The differential size distribution follows a clas- 
sical power law d n(a) oc a K da, from a minimum size to 
a maximum size a mm . The minimum grain size <7 m i n and the 
power law exponent k are free parameters that are constrained 
by the SED fitting (Sect. 4). The maximum grain size a ma x is 
a fixed parameter, sufficiently large so that it does not affect 
the light emission properties at the wavelengths considered (we 
chose flmax ~ 8 mm). For standard k values (close to -3.5 or 
smaller), the mass of solid material in a debris disk can be “arbi- 
trarily” large depending on the assumed higher cut-off value of 
the size distribution (that can extend up to planetesimal sizes). 
In this article, we define the dust mass A/d ust as the total mass of 


material contained in grains smaller than 1 mm. This definition 
has been followed within the GASPS consortium because it al- 
lows direct comparison between different studies independently 
front the model assumptions. 

We adopt a complex grain model that bears some similari- 
ties to those succesfully employed by several authors for fitting 
SEDs of young debris disks (Li & Greenberg 1998; Augereau 
et al. 1999; Li & Lunine 2003a,b), and for fitting scattered 
light disk colors (Kohler et al. 2008). We consider icy porous 
aggregates, with silicate cores and organic refractory mantles 
(or simplified versions of this model; see Sect. 4), to mimic 
fluffy particles made of essentially unalterred interstellar mate- 
rial (the so-called “cold-coagulation” dust model in Li & Lunine 
2003a, b). Assuming cosmic abundances, the volume of the car- 
bonaceous mantle Vc is estimated to be ~2.12 times that of the 
silicate core Vsi in the interstellar medium, a value similar to the 
one derived from the mass spectra of comet Halley in the solar 
system (Vc/Vsi ~ 2, Greenberg 1998). 

To calculate the optical properties, we use a number of repre- 
sentative materials with well established optical indexes. These 
are either derived from laboratory measurements or from obser- 
vational and analytical arguments. The materials are: (1) the as- 
tronomical silicates (astroSi) from Draine (2003), a theoretical 
material consisting of amorphous silicates polluted by carbona- 
ceous materials and other metals (p = 3.5 g cm -3 ); (2) the ACAR 
sample of amorphous carbon from Zubko et al. (1996) as organic 
refractories (p = 1.95 gem -3 ); (3) the H20-dominated ice from 
Li & Greenberg (1998) (amorphous ice, p = 1 .2 g cm -3 ); and (4) 
the crystalline water ice from Warren (1984) (p = 0.92 gent -3 ). 

The optical indexes of an aggregate are obtained using a 
standard effective medium theory (EMT), the Bruggeman mix- 
ing rule (Bohren & Huffman 1983), and following the method- 
ology of Li & Greenberg (1998) and Augereau et al. (1999). 
The optical properties are computed with Mie theory valid for 
hard spheres. The EMT-Mie combination has been shown to be 
a good approximation for the absorption and scattering efficien- 
cies, provided that grain shapes are close to ellipsoids such that 
the porosity can effectively be seen as the volume fraction of 
vacuum (Shen et al. 2009, and references therein). In fact, Shen 
et al. (2009) show that only the polarization and the phase func- 
tions are significantly sensitive to grain shapes. We define the 
composition of a grain through the volume fractions of its differ- 
ent components: 

Vsi Vc Vice _ Vvac 

Vq: = , Vr = , «ice = , r = (1) 

V mat v mat ’ V mat ’ V mat + V vac 

where Vs;, Vc, Vi ce and V vac are the total volume of astronomi- 
cal silicates, carbonaceous material, ice and vacuum in the grain, 
respectively, and V mat = Vsi + Vc + Vi ce is the total volume of 
solid material. V denotes the porosity. Should one of the materi- 
als within the grains reach its sublimation temperature, it is re- 
placed by vacuum and the optical properties are adjusted accord- 
ingly in the model. In most cases however, the disk is too cold for 
any material to sublimate. Figure 7 shows the mass opacity (ab- 
sorption cross section per unit mass) calculated using Eq. (4) of 
Draine (2006) for several prototypical dust compositions, with 
the size distributions given in Table 6. 

4. Modeling the dust continuum emission 

We use the Grater code (Augereau et al. 1999) to compute a grid 
of models with a goal of reproducing the HD 181327 spectral 
energy distribution. Grater is a SED, image and interferomet- 
ric data fitter, specially designed to efficiently model optically 


A17, page 6 of 16 


J. Lebreton et al. : An icy Kuiper belt around HD 181327 
Table 5. Parameter space explored with Grater , corresponding to a total of about 1 150 000 models. 



Description 

Parameter 

Explored range 

Values 

Distribution 


minimum grain size 

^min [/^m] 

0.01... 100 

77 

log 

Sect. 3.4 

size power law index 

K 

-6.0... -2.5 

50 

linear 

maximum grain size 

Umax [mm] 

7.8 

fixed 

- 


Disk mass up to 1 mm 

^dust [^©] 

>0 

fitted 

- 

1 AND 2-MATERIAL MODELS 

Sect. 4. 1.1 

Si+Amorphous C 

ACAR volume fraction Vq 

0.00... 1.00 

21 

linear 

Sect. 4.1.2 

Si+amorphous ice 

ice volume fraction t)j ce 

0.00 . . . 0.90 

19 

linear 

Sect. 4.1.2 

Si+crystalline ice 

ice volume fraction Cj ce 

0.00 . . . 0.90 

19 

linear 

Sect. 4.1.3 

Porous Si 

Porosity P 

0.00 . . . 0.99375“ 

24 

linear 

3 AND 4-MATERIAL MODELS 


Si 

astrosi volume fraction «si 

1/3 X (1 - Dice) 

- 

- 

Sect. 4.2 

+amorphous C 

am. carbon volume fraction Vc 

2/3 X (1 - U ice ) 

- 

- 

+amorphous ice 

ice volume fraction y ice 

0.00. . . 0.90 6 

12 

linear 


+ vacuum 

Porosity P 

0.00. . . 0.95“ 

18 

linear 


Notes. The volume fractions Vsi, Vc, v ; cei and the porosity P are defined in Sect. 3.4 (Eq. (1)). 19 linearly distributed points from 0.00 to 0.95 + ad- 
ditional high porosity points: 0.96875, 0.975, 0.9875, 0.99375. ih) linearly distributed points from 0.00 to 0.90 + additional points at 0.65 and 0.75. 
(c) 10 linearly distributed points from 0.00 to 0.90 + additional points at 0.45, 0.55, 0.625, 0.65, 0.675, 0.75, 0.85, 0.95. 
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Fig. 7. Mass opacity of spherical grains computed with the Mie theory 
and Bruggeman mixing rule for several prototypical compositions de- 
tailed in Table 6 and their associated best-fit size distributions. 


thin disks with parametric grain size and radial distributions, or 
distributions from dynamical simulations. The code accounts for 
both the scattered light and continuum emission of dust grains in 
thermal equilibrium with a star. 

In the present case, the spatial dust distribution is constrained 
by the resolved HST images (Sect. 3.1, Fig. 5). The surface den- 
sity profile derived from inverting the observed surface bright- 
ness was imposed in the model. Our study of the disk contin- 
uum emission focuses on the properties of the dust particles, 
namely their size distribution and composition. To that end, we 
computed large grids of models reflecting various possible dust 
properties, starting from the simplest composition (astronomi- 
cal silicates) and progressively increasing the model complex- 
ity until obtaining a good match to the SED. For each model, 
the dust mass Md us t and the quality of the fitting are obtained 
from linear least-square fit between the observed and the syn- 
thetic SEDs. The^ 2 listed in this section are obtained using the 
48 data points presented in Tables 2-4, and with the numbers 
of degrees of freedom listed in Table 6. We subsequently use a 



Fig. 8. Best fit models (Table 6) for several 1- or 2-material mixtures 
{top) and 3- to 4-material mixtures (bottom). Red crosses: photometric 
data (Fig. 6). Solid grey line: synthetic stellar spectrum. 


statistical (Bayesian) inference method to derive the probability 
distribution of each parameter (see Appendix A). 

In total we consider ~300 compositions, with 77 x 50 grain 
size distributions for each, leading to ~ 1 150 000 models. All the 
parameters of the grain models are summarized in Table 5. 
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Fig. 9. Bayesian probabilities for the models in Table 6: minimum grain size, size distribution power law index, and total dust mass. 


4.1. First steps toward a complex dust model 

We start by fitting the HD 181327 SED with a single- or two- 
material grain model, which show the need for a more complex 
dust model. We consider several mixtures of the materials pre- 
sented in Sect. 3.4 to assess the broad properties of the grains. 
The model then includes 3 or 4 free parameters: (1) the mini- 
mum grain size a m i„; (2) the slope of the size distribution k ; (3) 
the total mass of the disk Mdust; and in the case of 2-component 
grains: (4) the volume fraction of either carbon ( vc ) or ice (i>i ce ), 
or the porosity (P). The detailed results of the investigation are 
summarized in Table 6, the best-fitting SEDs are displayed in 
Fig. 8 (top panel) and probability curves are shown in Fig. 9. 

4.1 .1 . Silicates or organic refractories? 

The most basic grain model to consider is the case of com- 
pact spheres made of pure astronomical silicates (astroSi), a 
configuration that is often used by default in many model- 
ing approaches. We let a m j n , k and Md us t vary and we com- 
pute a grid of x 2 on which we apply a Bayesian analysis. 
We find that pure astronomical silicate grains would have 
r7 min = 0.99 ± 0.1 3 pra and k = -3.43 ± 0.04, leading to a total 
dust mass Mdust = 0.35 ± 0.03 M© for grains smaller than 1 mm 
(Table 6). These values are consistent with the ones typically 
found for debris disks, but they actually lead to very poor fits, as 
illustrated in Fig. 8 (top panel) where the model with the small- 
est^ 2 (Xr = 10.3) is displayed. Although the model is in fairly 
good agreement with the data over the full spectral range and 
agrees well with the sub-mm measurement, it fails to reproduce 
the shape of the emission peak at ~50-100/um and the mid- 
IR fluxes, shifting the dominant part of the disk to colder tem- 
peratures. It overestimates by 30-40% the excess at far-IR wave- 
lengths (A ~ 100-160 pm), in particular the HerschelfPACS 
measurements. 

We take a further step and consider mixtures of amorphous 
silicates and carbonaceous materials. We ran a Bayesian analysis 
and find that the actual best composition corresponds to pure 
carbon. Qualitatively, the effect on the SED is not significant (see 
Fig. 8, top panel), even if it leads to a better minimum chi-square 
(Xj = 7.75). This model gets closer to the measured fluxes at the 
peak emission (A ~ 50-100 pm), but does not provide an overall 
good fit. The addition of carbonaceous material slightly shifts 
the best fit parameters of the grain size distribution (Table 6), 
resulting in slightly larger a m ; n values and slightly steeper size 
distributions 

4.1.2. Icy grains? 

Chen et al. (2008) suggested that HD 181327 may contain a sig- 
nificant fraction of water in the form of crystalline ice. producing 


a broad peak in the combined Spitzer / IRS (5-35 pm) and MIPS- 
SED (60-77 pm) spectrum. To assess the presence of water ice 
in the debris disk of HD 181327, we model the disk with grains 
made of a mixture of astronomical silicates and water ice, either 
in the form of amorphous ice or crystalline ice. 

The results of the statistical analysis are documented in 
Table 6. The quality of the fit is only slightly improved com- 
pared to the pure silicate model: the best models have;y 2 = 8.90 
for amorphous ice, Xr = 10.4 for crystalline ice. In fact, the 
addition of water ice, crystalline or amorphous, does not mod- 
ify the general shape of the synthetic SED (see Fig. 8, top panel). 
The smallest ;p 2 calculated with amorphous ice is slightly smaller 
than the one of crystalline ice, but in both cases the values are 
too high to accept these models. 

4.1.3. Porous grains? 

As a final step in our preliminary study of the dust composi- 
tion, we consider a grid of models consisting of silicate grains 
with different porosities, from compact silicates to extremely 
fluffy silicate grains (P mm = 99.375%). The marginal proba- 
bility distribution for the porosity P reveals a bimodal proba- 
bility distribution: the most probable models are the ones with 
either no porosity (discussed in Sect. 4.1.1), or very high poros- 
ity (P = 97.5%). Both solutions are presented in Table 6. A 
visual inspection of Fig. 8 resolves the ambiguity. The com- 
pact grain model better reproduces the minimum and maximum 
wavelength ends of the SED, and in particular fits the sub-mm 
point, but it fails at reproducing the intermediate wavelength 
range around the emission peak. On the other hand, the highly 
porous grain model is orders of magnitude below the observed 
value at 870 pm, but it is within the error bars in the other regions 
of the SED, especially at Herschel wavelengths. In sum, none of 
these models lead to overall good fits to the whole SED. 

4.2. Construction of an advanced grain model 

From what precedes, we conclude that none of the 2-component 
mixtures considered are sufficient to provide a satisfactory fit to 
the SED. This is essentially a consequence of the disk surface 
density being known a priori in this model; it cannot be adjusted 
to compensate for the emission properties of the dust grains. We 
note, however, that the addition of one component to silicates, 
should it be carbonaceous material, amorphous ice, or poros- 
ity, provides some slight improvements. This suggests that their 
combined effect could produce a more realistic model. 

We evaluate the improvement provided by a 3-component 
mixture on the fitting of the SED. We showed in Sect. 4.1.1 that 
it is not clear how to disentangle the respective contribution of 
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Table 6. Best-fitting (Best, i.e. lowest^) and most probable (Proba, i.e. highest probability) parameters assuming different grain compositions 
(see text for the definition of the parameters). 


Composition 


Vsi 

«C 

^ice 

P 


K 

^dustCM®) 

Tr 

1 -MATERIAL (d.O.f. = 45) 









Si 

Proba : 

1.00 

- 

- 

- 

0.99 + 0.13 

-3.43 + 0.04 

0.35 ± 0.03 


Best : 

1.00 

- 

- 

- 

1.00 

-3.42 

0.35 

10.3 


Proba: 

- 

1.00 

- 

- 

1.62 + 0.03 

-3.65 ± 0.07 

0.15 + 0.03 



Best : 

- 

1.00 

- 

- 

1.62 

-3.64 

0.15 

7.75 

2-material (d.o.f. = 44) 









Si+C a) 

Proba: 

0.007f 

0.99 ± 0.04 

- 

- 

1.62 + 0.03 

-3.65 ± 0.07 

0.15 + 0.03 


Si+am.ice 

Proba: 

Best: 

o P 

- 

0.41 ± 0.05 
0.40 

- 

0.15 ± 0.13 
0.16 

-3.36 ± 0.07 
-3.35 

0.26 ± 0.05 
0.26 

8.90 

Si+cr.ice 

Proba: 

Best: 

0.67f 

0.65 

- 

0.33 ± 0.12 
0.35 

- 

0.29 ± 0.24 
0.26 

-3.35 ± 0.02 
-3.35 

0.27 ± 0.04 
0.27 

10.4 

Si+P 

Proba: 

0.94f 

- 

- 

0.057 ± 0.062 

0.98 ± 0.20 

-3.42 ± 0.05 

0.35 ± 0.04 


Best: 

0.025 

- 

- 

0.975 

0.010 

-3.28 

0.039 

11.9 

3-material (d.o.f. = 44) 









Si+C+P 

Proba: 

0.33t 

0.67t 

- 

0.95 ± 6xl0 -5 

9.45 ± 3.20 

-3.13 + 0.09 

0.042 ± 0.001 


Best: 

0.33f 

0.67f 

- 

0.95 

10.0 

-3.14 

0.042 

6.05 

Si+C+ice 

Proba: 

0.07f 

0.13f 

0.700 ± 0.002 

- 

0.78 ± 0.04 

-3.50 ± 0.03 

0.11 + 0.01 


Best: 

0.07f 

0.13t 

0.80 

- 

0.55 

-3.42 

0.11 

4.48 

4-material (d.o.f. = 43) 









Si+C+ice+P 

Proba: 

O.llf 

0.22f 

0.67 ± 0.07 

0.63 ±0.21 

0.81 ± 0.31 

-3.41 ± 0.09 

0.051 ± 0.016 


Best: 

0.12f 

0.23f 

0.65 

0.65 

0.89 

-3.42 

0.048 

1.64 


Notes. The most probable parameters are given with 3-cr uncertainties. The volume fractions marked with the f symbol are imposed knowing the 
other free volume fractions, “d.o.f.” means degrees of freedom (see Appendix A). (1) The best fit for the Si+C mixture is obtained for silicate-free 
grains. The best-fitting parameters for that solution are given in the 1 -material part of the table. 


silicates and organic materials by looking at the SED. Instead of 
including an additional free parameter, we consider grains (ap- 
proximately) in fixed Si/C ratio described in Sect. 3.4 ( uc/^si = 
2), and we restrict the discussion to the ice fraction and the 
porosity of the grains. We thus assume that grains are made of 
silicate and carbonaceous material, which are either icy (i>i ce pa- 
rameter) or porous (P parameter). The results are reported in 
Table 6. For porous ice-free grains, the best model is found for 
a porosity P = 95% and leads to some improvement to the fit 
with Xr = 6.05. If we now consider only compact icy grains, 
we can achieve a very reasonable^ = 4.48 with comparison to 
previous models, assuming an amorphous ice fraction Dj ce = 0.7. 
Overall, both the addition of either water ice or vacuum to sili- 
cate + carbonaceous grains do contribute to improving the fit to 
the SED of HD 181327. This can can be seen in Fig. 8 (bottom 
panel). Nevertheless, the fits are still unsatisfactory (x? much 
larger than 1), and justify an even more advanced dust model. 

We next consider a 4-component grain model inspired by 
the “cold coagulation dust model” by Li & Lunine (2003a,b), 
which has been successfully employed for modeling debris disks 
(Li & Greenberg 1998; Augereau et al. 1999; Li & Lunine 
2003a,b; Kohler et al. 2008). We keep the Si/C volume ratio 
fixed (vsi'.vc = 1:2) and let both the ice fraction and the porosity 
vary, in addition to the other usual free parameters in this study, 
namely the minimum grain size fl m in, the slope of the size dis- 
tribution k, and the total mass of the disk Mdust- This leads to 
a five-free-parameter disk model to fit the 48 data points in the 
SED. We compute three grids of models corresponding to the 
three surface density profiles shown in Fig. 5 for |gr H Gl = 0-30 
(i.e. the nominal and ±lcr profiles). The Bayesian analysis is 
then applied to the 3 grids separately and the probabilities for 
each point of the parameter space are calculated as the sum of 


the probabilities for the 3 profiles. The detailed results are docu- 
mented in Table 6. 

Figure 10 presents the results of the Bayesian analysis on the 
composition of the dust grains, with a 2D-map showing the prob- 
ability as a function of P and V[ ce for the grid models, after pro- 
jection of the probabilities over the other parameters. It clearly 
illustrates the need for a composition within a limited range of 
possible values. In particular, the fraction of H 2 0-dominated ices 
is well defined: the probability peak is obtained for an ice frac- 
tion t)j ce = 0.67±0.07. The porosity is not so strongly constrained 
but we find a peak of probability for P = 0.63 ± 0.21. The best 
fit model in terms of x 2 is found for grains made of 23% car- 
bonaceous materials, 12% silicates, 65% ice, and with a poros- 
ity P = 65%. As can be seen in the bottom panel of Fig. 8, for 
such grains, the synthetic SED very reliably describes the ob- 
served SED from the mid-IR to the sub-mm, with Xi = 1-64, a 
value significantly better than those obtained for simpler grain 
models. The grain size distribution derived from the statistical 
analysis has a m j n = 0.81 ± 0.31 fim and k = -3.41 ± 0.09, a 
power law index very consistent with that predicted for theoret- 
ical collisional cascades ( k = -3.5). The mass in grains smaller 
than 1 mm, Md us t = 0.051 ± 0.016A/ ffi , may look surprisingly 
small at first glance compared to the other models, but it results 
from the significant fractions of ice and vacuum in the grains, 
implying larger mass opacities (Fig. 7). The temperature of the 
grains at the surface density peak position (89.5 AU) ranges from 
86 K for the smallest grains to 40 K for ~8 mm-sized grains. 

To summarize, our approach, based on a fixed surface den- 
sity profile constrained by resolved imagery of the disk, clearly 
highlights the need to consider the presence of both ice and 
porosity to find an overall good fit to the SED. Icy grains are 
naturally expected at such a large distance from a solar-type 
star and the fluffiness of the grains is likely a by-product of 
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Porosity 


Fig. 10. Bayesian probabilities for the grain composition in the case of 4-material dust mixtures. 


their formation process. We applied a similar cold-coagulation 
dust model as Li & Lunine (2003a) did fitting the very young 
HD 14 1569 A debris disk, except that we used no a priori as- 
sumption on the ice fraction, but with arguing the grain com- 
position must reflect the primordial abundances of condensable 
elements, they chose Dice ~ 30% for cold regions and they find 
a best-fitting porosity P ~ 0.8. Hence we propose a dust model 
that is especially ice-rich with a degree of porosity that is not 
unexpected. Further comparisons with the literature and impli- 
cations on the nature of the dust are discussed in more details in 
Sect. 6. 


5. Limits on the gas content 

At 12 Myr-old, HD 181327 offers a chance to assess the amount 
of gas in a young debris disk, at a time when, according to planet 
formation theories, the accretion of gas by giant planets has al- 
ready occurred while the formation of rocky planets might be 
ongoing. 

We use the upper limits on the oxygen and ionized car- 
bon fine-structure lines ([O I] at 63 and 145 yum, and [C II] at 
158 /rm, see Table 3) to constrain the amount of gas present in 
the disk. We solve the photochemistry and energy balance in the 
disk using the photochemical code ProDiMo . The abundance 
of 7 1 species is computed in steady-state, under the influence of 
adsorption, desorption (thermal, photo and cosmic-ray induced), 
and molecular hydrogen formation on grain surfaces. 

The gas density profile is supposed to match the one of the 
dust, as presented in Sect. 3.1. The total mass of solid (defined 
as the mass in grains from 0.81 yum to 8 mm) is constrained by 
the best dust model of Sect. 4.2 (M so iid = 4.9 x 10 -7 M Q = 
0.164 M e ) and we let the total atomic+molecular gas mass M gas 
vary from a gas-poor disk (M gas /M so ua = 10 -2 ) to a gas-rich 
disk (M gas /Af so iid = 10 3 ). The parameters of the model are sum- 
marized in Table 7. One of the main parameters controlling the 
gas energy balance is the amount of PAHs (/pah), which remains 
unconstrained in the disk of HD 181327 because their abundance 
is too small to produce detectable PAH emission features in the 
IR. We adopted two “extreme” values for /pah'- a low PAH abun- 
dance (/pah = 10 -5 ) and a high PAH abundance (/pah = 0.1) 4 . 

The non-LTE line fluxes for the [OI], [CII], CO 7 = 3 
— > 2, and CO J = 2 — > 1 transitions are given in Fig. 1 1 and 

3 Woitke et al. (2009) provide a detailed description of the code, with 
subsequent improvements described in Kamp et al. (2010), Thi et al. 
(2011), and Woitke et al. (2011). 

4 The high PAH abundance is an upper-limit to the amount of PAHs in 
HerbigAe disks, where it is found to be of the order of 0.1-0.01 (e.g., 
Meeus et al. 2010, / PAH = 0.03 for HD 169142). 


Appendix B. Among all the lines, only the CO J = 3 — > 2 and 
J = 2 — > 1 line fluxes behave monotonically, increasing in flux 
with increasing disk mass. The reasons for the apparent erratic 
behaviour of the fine-structure lines come from the non-linear 
nature of several physico-chemical phenomena. For example, 
the [C II] flux depends on the abundance of C + and on the exci- 
tation of the first energy level. The C + abundance first increases 
with total gas mass until CO self-shielding against photodisso- 
ciation is efficient enough such that a significant fraction of C + 
is converted into neutral C and CO. The presence of CO also 
changes the heating and cooling balance, since the CO absorp- 
tion of IR stellar photons becomes the main source of heating. 
The change from a disk with no CO self- shielding to a CO self- 
shielded disk manifests itself as a jump in CO fluxes between a 
disk with gas-to-solid ratio of 1 to a disk with gas-to-solid ra- 
tio of 10. The abundance of PAHs has a weak influence on the 
gas temperature because the small UV flux from a F5/F6 star 
without excess accretion luminosity cannot ionize PAHs. On the 
other extreme, when the gas mass is very small (at gas-to-solid 
mass ratio below 1), the PAH and molecular absorption heating 
are nonexistent. Photoelectric heating by large grains dominates 
and explains why the line fluxes are similar for the high and low 
PAH abundance models. 

Line fluxes from models with /pah = 10 -5 are consistent 
with the PACS upper limits for all gas-to-solid mass ratios (up- 
per part of Table B.l). Only in the case / PA h = 0.1, the disk 
gas mass is constrained by [OI] at 63 yum to be below ~17 M ® 
(A/ gas /M S oiid 55 10 2 , lower part of Table B.l). 

Eventually, even if the integration of the [O I] 63 and 
145 /im, and [C II] lines we performed was not sufficient to de- 
rive strong constraints on the gas-to-solid mass ratio in the disk, 
the combination of the GRaTeR and ProDiMo codes proved effi- 
cient to make predictions on the fluxes in the CO J = 3 — > 2 and 
J = 2 — > 1 lines. These lines will be detectable for instance by 
ALMA, finally allowing to assess the gas content of the debris 
disk. 

6. Additional constraints on the disk 

The surface density profile used in this study comes from a new 
reduction of the scattered light images of Schneider et al. (2006) 
obtained with HST/NICMOS (see Sect. 3.1). This ensures by 
definition that the shape of any synthetic scattered light sur- 
face brightness profile in this study perfectly agrees with the 
observed profile at 1.1 /mi. However, the consistency is limited 
to the shape, and the agreement in terms of absolute flux is not 
guaranteed. This is because our fitting approach in Sect. 4 con- 
sists of adjusting the disk mass to get a best fit to the SED, which 
is equivalent to scaling up or down the surface density profile. A 
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Fig. 11. Expected gas line fluxes as a function of gas-to-solid mass ratio 
for either a low-PAH abundance (thin lines) or a high-PAH abundance 
(thick lines). The horizontal lines denote the PACS 3cr upper-limits and 
do not include the 30% flux calibration uncertainty. The gas-to-solid 
mass ratios assume a solid mass of 0.164 M a . 

Table 7. Setup of ProDiMo for the gas disk model. 


Parameter 

Value 

Elemental abundances 

ISM-like 

UV field 

1 Habing field" 

Cosmic-ray ionization rate ( 

1.7 x 10- 17 s -li 

Non-thermal line width dv 

0.15 km s -1 

/pah c 

io - 5 -io -1 


Notes. <o) Dalgamo (2006), (6) Standard interstellar UV field, (c) / PA h 
scales with the gas mass and is defined with respect to an interstellar 
abundance of 3 x 10 -7 . 


reason for adopting this strategy comes from the fact that the 
scattered light images do not directly provide the surface density 
E(r), but its product by the mean scattering cross section <x sca 
at the observing wavelength (1.1 pm). The fitting process to the 
SED, on the other hand, yields an independent, absolute scaling 
of the surface density profile, thus allowing one to obtain an ini- 
tial estimate of cr sca . We find cr sca = 3.7 pm 2 in the NICMOS 
1.1 pm band. 

A second estimate of cr sca is obtained by calculating 
the theoretical scattering cross section for compact spherical 
grains given the grain properties (composition, size distribution) 
that were inferred in Sect. 4: cr sca (d) = J_ m “ 7m 2 Q sca (a, A)dn(a) 
where Q sca is the scattering efficiency at wavelength A and for 
a grain radius a. For the best model (in the Bayesian sense, 
Table 6) we find <x sca = 16.8 pm 2 at 1 . 1 pm, a factor of about 4.5 
higher than the first estimate. This implies a peak flux density of 
3.14mJy arcsec -1 inconsistent with the 0.76 mjy arcsec -1 mea- 
sured in the NICMOS profile. In other words, this means that 
our model predicts an albedo (at = 59% at 1 .1 pm) that is a fac- 
tor of about 4.5 larger than the observed value. This suggests 
that smaller grains than those required to fit the SED dominate 
the scattering behavior of the disk. It is worth noting that such 
a conclusion is not specific to HD 181327, but is also observed 
in other debris disks. Recent examples include HD 207129 for 
which Krist et al. (2010) predict an albedo (based on dust proper- 
ties also obtained from a fit to the SED) that is an order of magni- 
tude larger than the value suggested by the scattered light image. 


A similar inconsistency is noticed for HD 92945 by Golimowski 
et al. (2011) who point out that different observing diagnostics 
suggest different minimum grain size. 

It is important to keep in mind that although our model sug- 
gests the grains are very porous and hence likely to have an 
irregular, fractal structure, we used hard spheres (Mie theory) 
to calculate both the thermal emission and scattered light opti- 
cal properties. Unlike the thermal behavior, light scattering will 
likely strongly depend on the exact shape of the dust aggregates, 
and on the size of their elemental constituents, two features that 
are not accounted for by the Mie theory (Voshchinnikov et al. 
2007). Another piece of evidence that Mie theory might be a 
limiting factor for a self-consistent model of HD 181327 comes 
from the observed scattering asymmetry parameter qhg that is 
a measure of the direction of the light scattered by a particle. 
For our best model (Sect. 4.2) the Mie theory yields a mean 
scattering asymmetry parameter I^hgI = 0.95 at A = 1 .1 pun, 
much larger than the |§ihgI = 0.30 measured on the NICMOS 
image. Only a population of much smaller particles could pro- 
duce such a low scattering asymmetry parameter if we stick to 
the Mie formalism, as already noted by Schneider et al. (2006). 
An analogous problem arises from the modeling of the debris 
disk of HD 92945 (Golimowski et al. 2011). Our micron-sized 
grains are highly anisotropic, essentially due to their spherical 
shapes. Incorporating vacuum homogeneously inside the grains 
does not suffice to decrease |gf H d and a more accurate description 
of the shape of fluffy aggregates is eventually needed. In fact, the 
scattering efficiencies and asymmetry parameters of fluffy aggre- 
gates are known to be correctly described by an EMT-Mie model 
(porous spheres) only in the Rayleigh regime, i.e. when the size 
of the inclusion (vacuum here) is small compared to the wave- 
length (see e.g. Voshchinnikov et al. 2007). This is not the case 
for the relatively large grain sizes and porosity when the wave- 
length becomes small. 

6.1. A narrower parent belt? 

Unlike scattering, we expect the continuum thermal emission to 
be largely independent of the exact shape of the grains. Our best 
model of HD 181327’s disk produces a sub-mm flux of 39 mjy 
at 870/rm, consistent within 2cr with the LABOCA measure- 
ment (52 ± 6mJy). The ATCA measurement at 3.2mm was not 
model-fitted, because of large uncertainties on the optical con- 
stants at millimetre wavelengths. As a sanity check, we calcu- 
late the flux in the ATCA band using a power-law extrapolation 
of the SED. In the far-IR to mm regime, the model emission de- 
clines as /U 2 ' 51±0 ' 165 . This results in a flux equal to l.b^ jmJy, 
which is within ~3 <x compatible with the measured flux of 
0.72 ± 0.25 mjy. Thus, the model remains valid for the longest 
wavelength observations (coldest emitting particles), showing 
that the ring resolved in the near-IR with the HST is compatible 
with the spatial distribution of the largest grains, at least at first 
order. 

A limit of our approach is the assumption that the grain size 
distribution is the same everywhere in the disk, while it is known 
that only the smallest grains can populate the outer disk due 
to radiation pressure. The maximum size and power law slope 
of the dust size distribution should vary with the distance from 
the star (e.g. Fig. 4 in Augereau et al. 2001), whereas the model 
keeps both large grains in the outer disk and a constant slope for 


5 Note that, at longer wavelengths, the emission seems to falls-off 
faster than a simple power-law. 
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the size distribution. Such a radial variation of the size distribu- 
tion is suggested by a comparison of a synthetic image at 70 jam 
with the PACS 70yum map (Sect. 2.1). The synthetic image ap- 
pears broader than the observed one with a FWHM of 7.68", 
versus 6.44" for the PACS image. Our interpretation is that the 
PACS band reveals larger grains that experience less the force 
due to the radiation pressure and thus are more bound to the par- 
ent belt. Assuming all grains are equally distributed results in 
putting too many large grains in the outer disk: a more complex 
model should account for the true dynamics of the particles. 

6.2. Emission from inside the belt? 

The region of the disk interior to 75 AU cannot be reliably as- 
sessed from the NICMOS data because of the coronagraph and 
artefacts from the PSF-substraction interior to 1.5" (78 AU). 
Therefore we cannot exclude the existence of an additional dust 
component interior to the main belt, that would be invisible in 
the surface density profile we use. If micron-sized grains pop- 
ulate the inner disk, they would radiate mainly in the mid-IR, 
which might affect the results of our study. From the Gemini/T- 
ReCS images of Chen et al. (2008), we estimate that less than 
~15% of the light emitted by the disk at 18.3 yum originates from 
a region interior to 1" (Fig. 2 in Chen et al. 2008). 

We perform a simple test of the robustness of our results. 
Suppose that a belt of black body dust particles is orbiting the 
star at a distance r < 1". We calculate a black body spectrum 
than peaks at 18.3 gm ( T = 150 K) and we adjust it such that 
its maximum flux represents 15% of the 53mJy attributed to 
the disk at that wavelength. Then we subtract the black body 
flux from the observed SED of Table 4 to simulate an obser- 
vation of the system excluding the hypothetical inner compo- 
nent. We reproduce the approach of Sect. 4.2 with this new SED. 
Interestingly the results of the Bayesian analysis are not signif- 
icantly impacted. To account for the deficit of flux at mid-IR 
wavelengths, only the minimum grain size needs to be slightly 
adjusted to eliminate the smallest grains: a m ; n = 1.0 yum and 
k = -3.43, while we have £tmin = 0.89 /rm and k = -3.42 for 
our best fit (Table 6). The composition remains unchanged for 
the best model (xj: = 1.52) and no change is induced that would 
improve the consistency of the model with the scattered light flux 
(Sect. 6). Consequently we cannot reject the possible existence 
of an inner component to the disk, but none of our conclusions 
would suffer from its existence. 


7. Grain dynamics and ice origin 


7.1 . Blowout size 


Grains in the HD 181327 debris disk for which radiation pres- 
sure overcomes gravity are put on unbound orbits and eventually 
ejected from the system on a short time scale (see for instance 
Krivov 2010, and ref. therein). The balance between radiation 
pressure (F P ) and stellar gravity ( Fq ) forces is evaluated with 
the distance independent, but size-dependent, /? pr factor, given 
by 


n _ l£p[ _ O- (6pr)E* 

~ |F g | ~~ m AncGM * 


( 2 ) 


where (<2 pr ) is the pressure efficiency averaged over the stellar 
spectrum, m the grain mass, <x the grain geometrical cross sec- 
tion, c the speed of light, G the gravitational constant, and L* 
and M* documented in Table 1 . Assuming that the grains are 



grain size a [gm] 

Fig. 12. yS pr = |Fp|/|Fol as a function of grain radius. The horizontal 
dashed line denotes the limit between bound (0 < 1/2) and unbound 
(0 > 1/2) grains. 


spherical, Eq. (2) yields a critical size below which the grains 
originating from the parent belt are ejected from the system: 
flbiow = aiflpr = 1/2) (assuming circular orbits for the parent bod- 
ies). 

In Fig. 12, we represent the J3 pr ratios as a function of grain 
radius for different compositions. Our best 4-component grain 
model has a blowout size «biow = 4.83 yum, that is 5 times larger 
than the inferred minimum grain size in the system. Models im- 
plying no porosity nor carbon, on the other hand, have blowout 
sizes of the order of the inferred value of a m j n , more consis- 
tent with the theoretical argument that grains below the blowout 
size should be expelled from the system. The response of sil- 
icate grains to the stellar irradiation and gravity fields is only 
slightly affected by the presence of ice, leading to similar shape 
and moderated values for y 8 pr (a). Organic refractory grains on 
the other hand are highly sensitive to the radiation pressure and 
they can reach high/? pr values. Figure 12 shows that porosity has 
the effect of flattening the y 6 pl (a) function, pushing the critical 
value y8 pr = 1 /2 to larger grains. Hence the addition of porosity 
to a mixture of carbon and silicates produces a strong effect: yS pr 
not only reaches a high value at its maximum due to carbon, but 
the yS pr (a) function decreases slowly such that the critical value 
yS pr = 1/2 is reached only for large grain sizes. 

Decreasing the amount of carbon in the grains (it was im- 
posed somehow arbitrarily, see Sect. 4.2) would mechanically 
decrease their blowout size. Pure silicate grains with the same 
porosity and ice content have flbiow = 0.89yum, but they do not 
provide good fits to the SED. We computed additional grid of 
models with alternative vc/vsi values. The results confirm the 
conclusion of Sect. 4.1.1: all Si+C+HjO+porosity models lead 
to similarly good fits to the SED, as long as a sufficient fraction 
of carbon is considered. For instance, if silicates represent only 
1 0% of the C+Si mixture, then Di ce remains at 0.65 and the poros- 
ity is reduced to 55%, with Mi ust = 0.043 M®, a m ; n = 0.89 /tm 
and Xr = 1 -66. All these values remain within the uncertainties 
of the original model. Models involving carbon-poor grains do 
not provide good fits to the SED, but as long as vc/vsi > 0.5, the 
results remain unchanged. In summary, the discrepancy between 
flmin and ttbiow does not lie on the hypothesis we made on the 
composition of the grains. 

Our estimation of flbiow is based on the assumption that 
the grains are spherical and homogeneous, but the actual re- 
sponse of inhomogeneous aggregates to radiation pressure can- 
not be that simple. Saija et al. (2003) for instance show that the 
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fluffiness of the grains tends to decrease their radiation pressure 
cross-section, as if aggregates behaved like their elemental con- 
stituents. This would make grains below the theoretical blowout 
size less sensitive to radiation pressure, and allow them to remain 
bound to the star. Thus the apparent inconsistency may be inter- 
preted as another signpost of the complex structure of the grains. 
We emphasize that fixing a m ; n as equal to «biow calculated with 
the (EMT-)Mie theory would not be a relevant approach when 
dealing with porous grains. 

Another possibility is that grains are not in dynamical equi- 
librium (i.e. an imbalance between production and ejection), fol- 
lowing for instance a major collisional event, such as an asteroid 
break-up. An excess of dust grains from all sizes would have 
been produced, including grains below the blowout size. These 
grains could have survived destructive collisions and radiation 
pressure blowing for long enough to be observed. Of course, ob- 
serving such an event is unlikely, because dust typical lifetime is 
only of the order of 1000 years, which is small with comparison 
to the age of the star (see Sect. 7.2). But this “rare event hypoth- 
esis” would be consistent the scenario proposed by Chen et al. 
(2008) to justify the brightness asymmetry they claim to observe 
in thermal light. 

7.2. Time scales and origin of the ice 

Collisions occurring between grains in a debris disk are a key 
process as they generally result in the production of smaller 
dust particles, either by fragmentation of the target, or exca- 
vation. Any phenomenon longer than the characteristic colli- 
sional time scale will be considered negligible in the evolution 
of the dust population. Based on the principle that any particle 
will cross the disk midplane twice during one orbit, Thebault 
& Augereau (2007) proposed a simplified relation for the colli- 
sional time scale, valid for the smallest grains: f co i = f 0 rbit/TT°, 
with f orbit the orbital period at the surface density peak, and -r[ eo 
the geometrical vertical optical depth. Here the orbital period 
at 89.5 AU is ~850years, leading to t co \ ~ 7600years for the 
best model. This is more than three orders of magnitude smaller 
than the age of the system and two orders of magnitude smaller 
than the characteristic time-scale for Poynting-Robertson drag: 
f pr = 0.1 Myr for /? pr = 1/2 grains (Wyatt 2008). The evolution 
of the HD 181327 debris disk is clearly dominated by collisions 
requiring a mechanism to generate dust on a 10 3 — 10 4 years char- 
acteristic time scale. 

Attempts to model time-dependent chemistry with ProDiMo 
reveals that the grain water ice mantle would be photo- 
evaporated in a few hundred years (for ~1 /rm grains), to a few 
thousand years (for ~10/rm grains), a similar (or shorter) time 
scale to the dust mutual collision, requiring as well an ice re- 
processing mechanism. This provides a consistent view that the 
icy dust grains are maintained through collisions in a reservoir 
of icy planetesimals. We highlight that, while detecting the mid- 
infrared spectral features in the emission of a debris disk is very 
challenging - unambiguous detection have only been made in 
brighter circumstellar disks (see, e.g. Malfait et al. 1999) - the 
sole analysis of the continuum SED alone can provide clues on 
the amount of ice in dust disks, thus on the nature of the parent 
bodies, often referred to as the leftovers of planet formation. 

In the solar system Kuiper belt, the study of icy planetesi- 
mals, Kuiper belt objects (KBOs) in particular, is a prolific topic 
(for a recent review see e.g. Jewitt 2010). The presence of ice, 
together with the porous nature of the objects, is inferred from 
estimates of their densities. Unambiguous detection of water ice 
has been made possible by the spectroscopic near-IR signature 


of crystalline ice, for instance by Delsanti et al. (2010) who 
completed their study of the plutino Orcus by radiative transfer 
considerations concluding that it contains a mixture dominated 
by amorphous and crystalline water ice. Confronted with space 
weathering, crystalline ice should be amorphidized in a short 
timescale. This implies that ice needs to be stored (and possibly 
produced) in the interior of the object before being brought to the 
surface, possibly through cryovolcanic events. The nature of wa- 
ter ice in KBOs remains unclear, but the vast majority of these 
objects are too small to experience such geological processes 
(Dumas et al. 2011, and references therein). Thus KBOs covered 
by crystalline ice can only constitute particular examples, and 
the water ice covering the collisionally produced dust in the solar 
system is more likely amorphous. The young age of HD 181327 
allows the possibility that primordial crystalline ice has survived 
in its debris disk. Our study tilts in favor of the amorphous nature 
of that ice, suggesting that either ice was predominantly incorpo- 
rated in an amorphous form in HD 181327’s Kuiper belt objects, 
or that the amorphization timescale is shorter than ~10 Myr in 
this distant belt. 

8. Discussion 

The properties of the HD 181327 debris disk can be compared 
to those of the disks surrounding two of its co-moving “sis- 
ter” stars, ft Pictoris (A3V, 1.7 M 0 ) and AU Microscopii (M1V, 
0.5 M 0 ), which formed in the same molecular cloud from the 
same material and are coeval. Of course comparisons are lim- 
ited by the specific models used to study each object. The rel- 
atively large porosity found here should be put in the context 
of the grain properties found for AU Mic (Graham et al. 2007; 
Fitzgerald et al. 2007). Assuming a Si+C+H20 mixture, only 
highly porous grains CP = 80%), reproduces the scattered light 
profiles (exhibiting strong forward scattering), the high polar- 
ization and the thermal SED. In that model, the Mie theory 
remains valid to represent the scattering properties of porous 
grains, because the grains filling the outer disk are very small 
and stay in the Rayleigh regime. A total dust mass of 0.01 M $ 
is found, that essentially reflects the contribution from the in- 
ner mm-sized grains, concentrated in a ~40 AU belt. A simi- 
lar porous grain model was used by Li & Greenberg (1998) to 
model /3 Pictoris yielding a dust mass ~0.4 A/®, contained in 
the ~100 AU-wide outer disk. Overall, “cold coagulation-like” 
dust models prove effective at reproducing observations of these 
young debris disks. A trend appears for increasing dust mass 
with spectral type as could be intuitively expected, but despite 
this difference of scale, the three stars are orbited by a similar 
material. 

Given the young age of these objects, we are investigating 
the conditions holding within disks just after the planetary em- 
bryos may have formed, the gas dispersed, and possible accre- 
tion onto outer giant planets happened. To explain the belt-like 
structure of the HD 181327 disk, the most straightforward expla- 
nation is to invoke the dynamical influence of a planet inside the 
belt. The Wisdom (1980) analytical criteria gives the width of 
the chaotic zone a planet mean-motion resonances would create. 
Assuming circular orbits for the planetesimals, as suggested by 
the sharpness of the belt, this yields a range of possible values for 
the planet semi-major axis depending on its mass. Assuming the 
inner edge of the disk is located at 78 AU from the star, the yet 
undiscovered planet would lie from 73 AU for a Neptune-mass 
planet, down to 59 AU for a 10 Jupiter-mass planet. 

The material at that stage appears as very hydrated provid- 
ing a considerable reservoir of water in the outer system. The 
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existence of such a reservoir is hypothesized in scenarios of the 
solar system formation to explain the unexpectedly large water 
content on Earth for instance (Morbidelli et al. 2000). Indeed, 
terrestrial planets at this age (12 Myr) are likely too hot to hold 
on to water and later impacts of icy comet-like bodies from 
the outer regions are a possible mechanism for delivering wa- 
ter and volatiles to planets in the habitable zone. Thus, not only 
HD 181327’s debris disk reveals that the system has already 
grown planetary embryos - including maybe the bricks of rocky 
planets in the inner disk - but it likely holds the conditions that 
were prevailing before the solar system late heavy bombardment 
that depleted the Kuiper belt and determined terrestrial planets 
history. 

9. Summary and conclusions 

We have presented new Herschel / PACS and ATCA imaging of 
the HD 181327 debris disk, a 12 Myr old star belonging to the 
p Pictoris moving group. The 70 and 100 yum images as well as 
the 3.2 mm observation reveal a marginally resolved extended 
emission. Fluxes in the three PACS photometric channels were 
derived and complemented by continuum fluxes obtained from 
PACS spectroscopic observations. Together, these data provide 
precise information on the continuum thermal emission of the 
debris disk, but detailed models such as the one developed in this 
study are required to interpret the precise nature of the emitting 
material. Degeneracies between the grain properties and disk’s 
spatial structure are broken in our model by making use of a 
reprocessed HST/NICMOS coronagraphic image revealing that 
the dust is bound to a narrow belt peaking at 89.5 AU from the 
central star. We focus our study on the dust properties (com- 
position, size distribution) using a radiative transfer code, and 
run a large grid of debris disk models for subsequent statistical 
(Bayesian) analysis. We conclude that the thermal light spectrum 
of the disk is well reproduced assuming a mixture of silicates, 
carbonaceous material and amorphous ice, as well as vacuum to 
mimic a significant porosity. The observed SED can be very re- 
liably reproduced by a population of grains larger than ~ 1 /a m, 
which follows a power-law distribution matching the ideal colli- 
sional equilibrium case and represents ~0.05 M e (1.5x 10 -7 M 0 ) 
in grains with sizes up to 1 mm. Extending the size distribution 
to the parent bodies of the observed dust grains leads to a mass 
of ~50 M e for planetesimals up to 1 km. The composition we 
infer for the dust grains must be a good representation of a con- 
siderable reservoir of material hidden in unseen planetesimals. 

In this study, we attempted to use HerschelfPACS spec- 
troscopy to detect for the first time far-infrared emission lines 
in the debris disk. The non-detections we obtain were used to 
investigate the gas content of the debris disk with our photo- 
chemical disk model ProDiMo based on the dust model found 
with GRaTeR. The main uncertainty lies in the amount of PAHs 
in the disk. The current non-detections of [OI] and [CII] lines 
alone do not provide unambiguous upper limits to the gas con- 
tent. Only in the case of a high PAH abundance, we can set an 
upper limit to the gas mass of ~17 M e (one Neptune mass). 
Coupling with other tracers, in particular with the CO lines ac- 
cessible to ALMA, offers much better prospects to reach lower 
limits on the low gas content of debris disks. However, steady- 
state chemistry is likely not ideally suited to debris disk model- 
ing, and in future studies more appropriate prescriptions, time- 
dependent chemistry in particular, will be explored. 

We note that both the low albedo estimated from the scat- 
tered light images and the dust dynamics are not well modelled 
by spherical grains filled with vacuum, and that the grains are 


likely complex fluffy aggregates. Though their thermal behavior 
is correctly described by the EMT-Mie theory, the scattering of 
light at short wavelengths and their sensitivity to radiation pres- 
sure require a more accurate description of the exact shape of 
the aggregates. Our statistical analysis highlights that ignoring 
the porosity of the grains or setting the minimum grain size as 
equal to the blowout size would result in incorrect interpreta- 
tion of the observations, in particular because the grain opacity 
is strongly dependent on the exact grain composition. 

We stress some necessary limitations of the models. First, we 
use a grain size distribution that is independent of the distance 
to the star, whereas it is known that only the smallest grains can 
populate the outer parts of the disk under the effect of radiation 
pressure. In this paper we discussed clues that the observed scat- 
tered light profile may not be a perfect representation of that seen 
in thermal light, by Herschel or ATCA. The long wavelength 
emission emanating from the largest grains might originate from 
an even narrower belt. Further studies will need to account for 
the radial differentiation of grain sizes with the objective to fit al- 
together resolved images from the near-IR to the millimetre do- 
main. Secondly, the Dohnanyi-like power-law size distribution 
is known to be an insufficient description; if the real distribu- 
tion implies wavy patterns as proposed by Thebault & Augereau 
(2007), then adjusting a power-law likely puts too many/too few 
small grains in the disk. Dynamical models are required to un- 
derstand better the general outcome of the collisional evolution 
of the disk. 
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Appendix A: Statistical analysis of the results 

The strength of the GRaTeR code, in comparison with an opti- 
mization algorithm that would provide a unique best-fit solution 
and possibly miss local minima, is that all the models are stored 
enabling fast and easy post-processing of the results. 

To find the best solution in the grid of models, we use a clas- 
sical least-square minimization approach with «f ree = 3 to 4 free 
parameters depending on the grain composition and N = 48, 
which results in 43 to 45 degrees of freedom (d.o.f.) (Table 6). 
The free parameters are k, a m i n and Md ust , and depending on the 
situation uc, V[ C t or P. Md us t is handled separately because it is 
automatically adjusted to reproduce the surface density profile 
in the code. 

Let a be the vector of parameters. Here a = {k.ci^.v^P). 
For each set of parameters, we can define the quantity: 


N 


x\a) = ^ 

i=l 


(D, - Mi ) 2 


(A.l) 


where D, and M, are the observed and modeled values, respec- 
tively, and a-j is the measured uncertainty. For each^- 2 value, we 
can define the likelihood of the data, i.e. the probability of the 
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Table B.l. Gas line fluxes. 


Gas-to-Solid 

[O I] 63 fim 

[O 11 145 A/m 

rc II] 158 yum 

CO J = 3 2 

CO J = 2 -> 1 

mass ratio 

(W m" 2 ) 

(W nr 2 ) 

(W m" 2 ) 

(W nr 2 ) 

(W nr 2 ) 



Herschel - PACS (3cr) 




9.8 x 10“ 18 

8.5 x 10“ 18 

7.0 x 10“ 18 

- 

- 



ProDiMo predictions with /pah = 10 5 


1000 

2.9 x 10“ 18 

1.6 x 10“ 18 

2.6 x lO" 20 

1.7 x 10“ 19 

5.6 x lO" 20 

100 

6.7 x 10“ 19 

1.2 x lO" 20 

4.1 x 10“ 20 

1.5 x 10“ 19 

5.3 x lO" 20 

10 

2.1 x 10“ 18 

8.2 x 10“ 21 

2.2 x 10“ 19 

1.3 x 10“ 19 

3.8 x lO" 20 

1 

2.4 x 10“ 18 

4.0 x lO" 20 

6.2 x 10“ 19 

3.4 x 10 -22 

6.1 x 10“ 23 

0.1 

6.4 x 10“ 19 

5.0 x lO" 20 

1.9 x 10“ 19 

1.0 x lO" 23 

1.6 x 10“ 24 

0.01 

2.8 x lO" 20 

2.9 x 10“ 21 

1.5 x lO" 20 

1.1 x 10“ 24 

2.7 x 10- 25 



ProDiMo predictions with /pah = 0. 1 


1000 

2.7 x 10- 17 

6.4 x 10“ 18 

4.1 x lO" 20 

4.8 x 10“ 19 

1.6 x 10“ 19 

100 

7.1 x 10“ 18 

2.0 x 10“ 19 

6.0 x lO" 20 

3.1 x 10“ 19 

1.0 x 10“ 19 

10 

3.2 x 10“ 18 

1.5 x lO" 20 

1.4 x 10“ 19 

1.4 x 10“ 19 

3.9 x lO" 20 

1 

2.4 x 10“ 18 

4.0 x lO" 20 

4.2 x 10“ 19 

2.6 x 10 -22 

4.6 x lO" 23 

0.1 

6.3 x 10“ 19 

4.8 x lO" 20 

1.8 x 10“ 19 

1.0 x 10 -23 

1.5 x 10“ 24 

0.01 

2.8 x lO" 20 

2.9 x 10“ 21 

1.5 x lO" 20 

1.1 x 10“ 24 

2.7 x 10 -25 


Notes. The gas to solid mass ratios assume a solid mass of 0.164 M e . The PACS 3-c upper limits do not include the 30% flux calibration 
uncertainty. 


data given the parameters of a model. This assumes the errors 
are distributed normally and are uncorrelated. 

P(D|a) = iexp(^^) (A.2) 

where C is a normalization constant and D = (D(i)}i=i..jv. In or- 
der to evaluate the effect of each individual parameter a t on the 
quality of the fit, we implement a Bayesian inference method. 
We define the posterior probability of a as the probability distri- 
bution of the parameters a given the data D. It has to be noted 
that this formulation implicitly assumes that the model is suited 
to describe the observation. Bayes’ theorem says that this quan- 
tity can be written: 


P(a\D) = 


P(a)P(D\a) 
P(D ) 


(A.3) 


where P{D) is the prior probability of the data under all possi- 
ble hypothesis (a normalization constant in practice), and P(a) 
is the prior probability distribution of the parameters a. Here we 
use uniform priors for all the parameters, meaning that we do 
not favor any realization of the parameters prior to obtaining the 
data. That fundamental theorem provides a simple way to inde- 
pendently assess the parameters of a model. Indeed, the marginal 
probability of each parameter a; can be computed regardless of 
the values of the other parameters aj t- i: 

P(a,\D) = f -f f-f P(a^\D)da 0 ...da^. (A.4) 

a 0 a i-l a i + 1 fl f'max 

Then we obtain a probability law as a function of the realization 
of each individual parameter. Finally, we estimate the param- 
eters by calculating the mathematical expectation and variance. 
These are the values presented in Table 6 where the uncertainties 
are3cr, 


E(iij) 


a ‘max 

J ' a,P(aj\D) 


(A.5) 


o- 2 (a,) = E(aj) - E(aif. 


(A. 6) 


Appendix B: Gas line fluxes (see Table B.l.) 
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